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1.0  INTRODUCTION 


The  Army  is  currently  developing  liquid  payloads  for  spin-stabilized 
projectiles.  A  liquid  filler  can  destabilize  the  flight  of  an  aeroballis- 
tically  well  designed  shell.  The  instability  mechanism  is  a  resonance 
between  the  coning  frequency  of  the  projectile  and  a  natural  vibration 
frequency  of  the  spinning  liquid.  Only  simplified  models  have  been  used 
in  the  past  to  study  this  liquid-induced,  projectile  instability.  These 
include  linear  models  and  axisymmetric  (independent  of  circumferential  or 
azimuthal  position  within  the  fluid)  finite-difference  solutions  to  the 
governing  equations,  the  Navi er-Stokes  equations  [Refs.  1,  2,  3] .  The 
linear  models  consist  of  closed  form  analytic  solutions  and/or  numerical 
solutions.  The  primary  deficiency  with  the  available  models  is  that  the 
projectile  yaw  must  be  quite  small.  In  fact,  experiments  to  verify  these 
modes  indicate  serious  nonlinear  effects  at  yaw  amplitudes  well  below  those 
common  to  projectiles  [Refs.  *♦,  5].  Hence,  large  yaw  effects  must  be 
incorporated  into  the  modeling  capabilities.  Also,  many  of  the  linear 
models  do  not  properly  treat  the  boundary  conditions  at  the  liquid/solid 
interfaces.  Corrections  are  required  within  these  models,  but  the  correc¬ 
tions  are  valid  only  for  high  Reynolds  numbers.  The  axisymmetric  Navier- 
Stokes  code  cannot  model  the  three-dimensional  disturbances  of  the  liquid. 
The  natural  frequencies  of  the  liquid  are  truly  three-dimensional  oscilla¬ 
tions  and  most  flow  problems  cannot  be  treated  by  an  axisymmetric  code.  A 
full  three-dimensional  solution  to  the  Nav ier-Stokes  equations  is  required. 
Viscous  effects  must  be  maintained  properly  in  all  coordinate  directions  to 
properly  model  the  rotating  fluid.  The  code  must  also  have  a  high  Reynolds 
number  capability.  This  requires  the  use  of  a  computational  algorithm  which 


is  more  complex  than  a  standard  finite-difference  solution. 

A  numerical  simulation  capability  developed  specifically  for  the  analy 
sis  of  fluid  flow  in  liquid-filled  projectiles  and  satisfying  the  above 
requirements  is  described  in  Ref.  6.  In  that  effort,  the  unsteady,  incomp¬ 
ressible  Navier-Stokes  equations  for  laminar  flow  are  solved  for  the  primi¬ 
tive  variables  without  recourse  to  linearization  or  simplification  of  the 
equations  of  motion.  The  finite  difference  approximations  to  the  govern¬ 
ing  equations  are  by  choice  either  first  or  second-order  time  accurate  and 
are  second-order  accurate  in  the  axial,  radial  and  azimuthal  directions 
(with  an  option  for  fourth-order  accuracy  in  the  azimuthal  direction).  The 
method  allows  imposition  of  arbitrary  body  motions  including  spin  and  pre¬ 
cession  and  the  corresponding  boundary  conditions  are  easily  and  directly 
prescribed.  The  finite  difference  equations  are  solved  using  an  implicit, 
approximate  factorization  procedure  that  permits  the  choice  of  reasonably 
large  time  steps  and  avoids  limitations  based  on  the  magnitude  of  Reynolds 
number.  Thus,  the  numerical  simulation  methodology  considered  provides  a 
complete,  and  flexible  framework  for  the  computational  analysis  of  fluid 
behavior  in  liquid  filled  projectiles. 

This  report  describes  progress  made  in  the  computattional  fluid  dyna¬ 
mics  of  fluid  filled  spinning  shells  beyond  the  capability  developed  in 
Ref.  6.  One  drawback  of  approximately  factorized  methods  is  the  error  due 
to  approximate  factorization  and  how  it  affects  the  accuracy  of  the  numer¬ 
ical  solution.  This  error  has  been  eliminated  in  the  current  effort. 
Rotating  frames  of  reference  have  been  utilized  resulting  in  simpler  codes 
for  simple  coning,  etc.  Simple  coning  motion  has  been  studied  along  with 
the  case  of  a  spinning  container  with  a  precessing  lid.  A  new  code  using 


the  axi symmetric  stream  function  formulation  has  also  been  developed. 

In  the  next  section,  the  derivation  of  implicit  upwind  schemes  without 
approximate  factorization  errors  is  outlined.  This  methodology  has  been  used 
in  all  the  codes  developed  in  the  current  effort.  The  third  section  deals 
with  the  development  of  the  axisymmetric  code  based  on  the  stream  function 
formulation.  The  development  of  unfactorized  implicit  schemes  for  the 
primitive  variable  formulation  is  outlined  in  the  fourth  section.  The 
fifth  section  deals  with  the  governing  equations  in  a  rotating  frame  of 
reference.  The  sixth  section  covers  the  application  of  the  method  to 
simple  coning  and  the  seventh  section  describes  the  application  to  a 
spinning  shell  with  a  precessing  lid.  Concluding  remarks  are  offered  in  the 
next  section  and  a  list  of  references  rounds  up  the  report. 
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2.0  IMPLICIT  UPWIND  SCHEMES  WITHOUT  APPROXIMATE  FACTORIZATION 


The  unsteady  Navier-Stokes  equations  can  be  divided  into  two  sets  of 
terms  accounting  in  turn  for  the  inviscid  and  viscous  behaviour  of  the  fluid. 
Along  with  the  time-derivative  terms,  the  inviscid  terms  constitute  a  hyper¬ 
bolic  system  of  equations  dominated  by  wave  phenomena  described  by  the  theory 
of  characteristics.  Implicit  finite  difference  schemes  for  the  unsteady 
Euler  and  Navier-Stokes  equations  have  for  the  most  part  used  central  space 
differencing  and  have  relied  upon  the  techniques  of  approximate  factorization 
or  fractional  steps  to  handle  multidimensions  [Refs.  7,  8],  Even  researchers 
using  various  upwind  schemes  (split-flux  method  or  Harten's  scheme)  [Refs.  9, 
10],  have  resorted  only  to  these  conventional  approximate  ways  of  splitting 
the  multidimensional  operators  into  one-d imens iona 1  operators  (which  are  then 
solved  efficiently  using  block-tridiagonal  elimination).  A  plus-minus 
splitting  scheme  [Ref.  ll]  has  also  been  attempted  for  the  split-flux  scheme 
leading  to  a  different  but  once  again  approximately  factored  implicit  scheme. 
In  contrast  to  the  above,  an  implicit  upwind  algorithm  is  used  in  the  current 
effort.  This  method  is  devoid  of  errors  of  any  kind  of  approximate  factori¬ 
zation.  The  advantages  of  using  such  an  algorithm  for  the  Navier-Stokes 
equations  for  incompressible  flow  will  be  explained  in  later  sections.  A 
brief  outline  of  the  fundamentals  underlying  the  new  approach  is  presented 
in  this  section.  For  a  more  complete  description  along  with  application  to 
the  compressible  Euler  equations,  Ref.  12  may  be  consulted. 

2.1  Theoretical  Framework 


It  will  be  shown  here  that  the  new  algorithm  owes  its  existence  to  the 
beneficial  properties  of  upwind  schemes  for  hyperbolic  systems  of  equations. 


It  cannot  be  constructed  for  central  differencing  schemes.  (These  remarks 
pertain  to  the  hyperbolic  part  of  the  Navier-Stokes  equations  only.  The 
viscous  terms  are  approximated  only  with  central  difference  formulae.)  At 
the  crux  of  the  new  method  is  the  observation  that  upwind  schemes  result 
in  a  diagonally  dominant  system  of  finite  difference  equations  governing 
the  change  of  dependent  variables  between  two  marching  (time)  steps.  Such 
a  system  can  then  be  solved  without  factorization  by  using  relaxation  iter¬ 
ations  between  the  two  time  levels.  The  Gauss-Seidel  iteration  technique 
lends  itself  naturally  in  this  context  as  will  soon  be  evident.  Any  other 
iteration  technique  will  work  if  it  is  applicable  to  a  diagonally  dominant 
system  of  equations  (thus  pointwise  Jacobi  iterations  will  also  suffice). 
The  choice  of  iteration  scheme  will  also  be  guided  by  the  internal  archi¬ 
tecture  of  the  computer  used  (parallel,  pipeline,  sequential  processors, 
etc.)  in  as  much  as  certain  relaxation  algoritms  lend  themselves  better  to 
certain  types  of  processing  to  achieve  faster  computational  speed. 

2.2  Implicit  Schemes 

We  are  concerned  here  with  hyperbolic  systems  of  equations  (including 
the  unsteady  Euler  equations)  in  mul t i -d imens ions .  With  the  addition  of 
centrally  differenced  viscous  terms,  the  methodology  is  easily  extended 
to  the  Navier-Stokes  equations.  Let  the  space  variables  be  x  and  y  and  the 
time  variable  t.  Only  two  spatial  dimensions  are  considered  here  but  all 
discussion  extends  in  a  straight  forward  manner  to  three  dimensions  (or  to 
one,  although  this  is  not  so  interesting).  We  will  be  considering  a 
quasi-1 inear  system  of  equations  of  the  type 
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qt  +  A(q)qx  +  B(q)q 


=  0 


(2.1) 


where  q  is  an  m-vector  and  A,  B  are  m  by  m  matrices.  For  simple 
presentation  of  the  underlying  ideas,  we  will  also  be  considering  in  more 
detail  the  linear  wave  equation  given  by 


u  +  au  +bu  =  0 

t  x  y 


(2.2) 


There  can  be  several  representations  and  versions  of  implicit  schemes, 
A  simple  implicit  scheme  may  be  constructed  for  the  equations  given  above 
by  employing  a  backward  discretization  of  the  time  derivative.  Thus,  the 
time  discretized  version  of  Eq.  2.1  is 


n+1  n 

q  -q 

At 


♦  An+,(q)q^' 


♦  Bn+,(q)qn+1 

y 


=  o 


(2.3) 


The  two  time  levels  have  been  indicated  by  n  and  n+1  .  Except  for  the 

case  of  linear  systems,  these  discretizations  obviously  result  in  nonlinear 

n+1 


equations  to  solve  for  q 
soon . 


The  spatial  discretizations  will  be  discussed 
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k- 

k- 

k- 
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By  linearizing  in  some  manner  the  Eqs.  2.3,  linear  implicit  prediction 
equations  can  be  constructed  for  qn+'  .  For  example,  we  can  consider 


n+1  n 

q  -q 

At 


.  An  n+I 
+  A  q 
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+  Bnqn+I  =  0 
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(2. A) 
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2.3  A  Simple  Upwind  Scheme 


Let  us  discuss  the  two-dimensional  wave  equation.  We  shall  see  how  a 
simple  upwind  scheme  for  it  results  in  a  diagonally  dominant  system  to  solve 
Then  we  shall  discuss  the  appropriateness  of  the  Gauss-Seidel  relaxation 
procedure . 

Assuming  a  and  b  are  positive,  backward  spatial  derivatives  can  be 

used  for  u  and  u  to  result  in  the  upwind  scheme 

x  y 


n+1  n 
u .  .  -u .  . 
J,k  J,k 
At 


+  a 


n+1  n+1 

u .  ,  -u .  ,  . 
J,k  J-i,k 

Ax 


+  b 


n+1  n+1 

J,k  J,k-1 
Ay 


0 


(2.5) 


This  scheme  is  first-order  accurate  in  space  and  time.  By  virtue  of  the 
linearity  of  the  governing  equation  (and  the  discretization),  Eq.  2.5  is 
a  linear  prediction  equation  for  un+'  . 


2.4  Diagonal  Dominance 

Rewriting  Eq.  2.5,  we  obtain 

t  1  a  _b\  n+1  _a_  n+1  _b  n+I  _  _l_  n 

At  Ax  +  Ay  Uj,k  Ax  U j - 1 , k  Ay  Uj,k.-1  At  Uj,k 

for  j  =  1 ,  .  . .  ,  J 

k  =  1 ,  ...,  K 


(2.6) 


It  is  clear  that  if  this  is  cast  in  matrix  notation  and  boundary  values  of 


uJ'J  and  uj+o  are  known»  t*ie  system  is  diagonally  dominant  by  rows  and 


columns.  In  other  words,  if  the  elements  of  the  matrix  are  d. 


(2.7a) 


Id.  ;|  >  E  I  d  j  J  for  j=l . J 

J,J  kflj  J,k 


>  I  |  d 


for  k=l , . . . ,K 


(2.7b) 


For  the  benefit  of  a  beginner  and  for  use  later  on  in  this  report,  the 
full  matrix  is  laid  out  in  Table  1  for  J=3  and  K=3.  It  is  clear  that  for 
large  J  and  K  ,  it  would  be  computationally  very  expensive  to  solve  the 
matrix  equation  by  direct  elimination.  (However,  in  the  special  case  of 
upwind  differencing  we  are  considering  for  the  case  a  and  b  positive, 
the  matrix  is  lower  triangular  and  can  be  efficiently  solved  by  direct  eli¬ 
mination;  the  above  remark  about  computational  expense  is  a  general  remark 
for  non-lower-triangular  and  non-upper-triangular  systems  of  difference 
equations).  Thus,  iteration  schemes  may  be  sought  for  to  solve  the  system 
of  equations.  Because  of  diagonal  dominance,  either  point-Jacobi  or 
Gauss-Seidel  iteration  procedures  may  be  used  to  iteratively  compute 
u.nt]  along  with  many  other  methods  valid  for  diagonally  dominant  systems. 


2.5  Jacobi  iteration 

Let  v!  ,  denote  successive  approximations  to  u?+'  .  The  Jacobi 

J  > K  J  » k 

(also  known  as  pointwise  or  simultaneous)  iteration  scheme  for  Eq.  2.6  may 
be  wr i tten  as 


lOiiu 

P  At 


i-1 

'  j  - 1  ,  k 


i-l  V 


(2.8a) 
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Matrix  form  of  implicit  upwind  scheme  for  wove 


q/(4x)  ,  r  =  b/(ay) 


While  convergence  is  guaranteed,  it  can  be  very  slow  for  pointwise  iteration. 


2.6  Gauss-Seide!  Iteration 


The  Gauss-Seidel  procedure  may  be  written  as 


1 

Vj,k 


,  u  •  i 

I  u A  + 

p  At 


q  vj-i  .k  +  r  v] ,k-i  > 


(2.9) 


It  is  assumed  here  that  the  notation  of  Eq.  2.9  implies  that  v^.  j  ^  and 
v!  ,  .  will  be  computed  before  v'.  .  and  thus  Eq.  2.9  is  an  explicit 

j , K- I  J 

formula  for  v!  .  .  Note  also  that  vi  ,  and  v!  n  are  known  from 

j,k  0,k  j,0 

boundary  conditions.  It  is  thus  clear  that  in  one  full  sweep  (or  fell  swoop) 

c  ,,  .  .  .  .  n+1  .  .  /.  first  iteration  n+1  v 

of  all  the  grid  points,  u.  ,  is  known  (i.e.  v.  ,  =  u.  .  ). 

J  J  J  »k 

To  make  the  procedure  even  clearer,  the  calculations  may  proceed  in  the 
following  order  of  grid  points  (corresponding  to  the  example  grid  of 
Table  1); 

1,1  to  2,1  to  3,1  to  or  simultaneously  with  1,2  to  1,3 

to  or  simultaneously  with  3,2 
to  3,3  .  (The  "simultaneously  with"  operations  can  be  exploited  for 

parallel  processing).  It  is  thus  clear  that  the  Gauss-Seidel  method  is 
actually  more  natural  for  the  wave  equation  than  for  a  linear  elliptic 
equat ion. 


to 

2,1 

to 

to 

2,2 

to 

2,3 

to 

3,3 
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process 

2.7  Notes  on  D i rec t i on  of  Sweep 


However,  it  is  also  clear  that  the  direction  of  sweep  is  very 
important  to  this  one  step  convergence.  Other  directions  of  sweep  implied 
for  example  by 
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2. 


J 


(backward  k  sweep  ,  forward  j  sweep) 


i  1  /  j,k  A  i-1  A  i  v 

j,k  p  At  M  j  -1  ,k  j.k-1' 

(backward  j  sweep  ,  forward  k  sweep) 


(2.11) 


will  lead  to  slow  convergence.  The  preferred  Gauss-Seidel  iteration,  Eq.  2.9 
is  actually  identical  to  a  direct  elimination  method  of  solution  of  the  mat¬ 
rix  system  of  equations  recognizing  the  fact  that  the  matrix  is  of  lower 
diagonal  form  and  very  sparse. 


2.8  Notes  on  Signal  Propagation  Directions 


The  constants  a  and  b  of  the  scalar  wave  equation  (Eq.  2.2)  signify 
signal  propagation  directions  by  their  signs  and  signal  propagation  speeds 
by  their  magnitudes.  When  a  and  b  are  both  positive,  signals  travel  to 
the  right  (along  the  positive  coordinate  direction)  in  both  x  and  y.  If 
a  or  b  had  been  negative,  it  would  have  been  appropriate  to  use  forward 
differences  for  the  space  der i vat i ve (s)  which  the  negative  constant(s) 
multiplied.  In  that  case,  the  directions  of  sweep  in  the  Gauss-Seidel  proce¬ 
dure  must  be  changed.  The  boundary  conditions  must  also  be  specified  along 
different  boundaries  than  those  natural  for  the  case  when  a  and  b  are 
positive. 


2.9  The  Variable  Coefficient  Linear  Wave  Equation 


Thus  far,  we  have  considered  a  and  b  to  be  of  the  same  sign  through¬ 
out  the  field.  In  such  situations,  one  of  the  four  different  combinations 


*■  w  -  -  »  -  • 


of  forward  or  backward  sweeps  in  j  or  k  will  result  in  one  step  conver¬ 
gence  of  the  Gauss-Seidel  procedure.  Let  us  consider  a  and  b  to  be  not 
constant,  but 


a(x,y,t)  ,  b  =  b (x , y, t) 


(2.12) 


A  generalized  upwind  scheme  for  Eq.  2.2  with  variable  coefficients  may  be 
written  in  such  a  manner  that  the  signs  of  a  and  b  are  automatically 
monitored  and  corresponding  one-sided  derivatives  assigned.  One  represent¬ 
ation  can  be 


.  n+1  n+1  n+1  n+1 

n+1  n  ii  u .  , -u .  .  .  ii  u .  i -u .  .  , 

u  -u  .  ;a+|aK  ,  j,k  j-1,k^  _  /a-]aK  ,  j  ,k  j+1  ,k^ 

At  ^2  M  Ax  ’  K  2  M  25  ' 

n+1  n+1  n+l_  n+l 

/ b+ 1  b j  %  ,uj  ,k~uj  ,k-K  _  /b-jbL  uj  ,k'uj,k+l  _ 
V  2  M  Ay  )  ^  2  '  Ay 


2.10  Cycling  Forward  and  Backward  Sweeps 


(2.13) 


The  following  sequence  of  two  Gauss-Seidel  (forward  and  backward) 
sweeps  is  expected  to  be  efficient  for  this  general  case. 


(_L  +  kL  +  M)v5 

'At  Ax  Ay  '  j.k 


(2.1*0 


/a+laK  il  /a-  a  \  i2  .  /b+lbL  il 

2  Ax  V j - 1 , k  "  2  Ax  )vj+l ,k  (2  Ay  )vj,k-l 


(b1kL)v>2  , 

'2  Ay  ’  j,k+l 


il  =  i  and  i2  =  i-1  for  odd  i,  il  =  i-1  and  i2  =  i  for  even  i 


When  a  and  b  are  of  the  same  sign  (+  or  -)  throughout  the  field,  one 
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cycle  (a  forward  and  a  backward  sweep)  is  enough  for  convergence.  In  the 
general  case,  when  a  and/or  b  change  sign  in  the  field,  one  cycle  will 
not  be  enough,  but  this  procedure  leads  to  a  quite  rapid  convergence. 


2.11  Second-Order  Accuracy 

We  have  so  far  discussed  only  algorithms  that  are  first-order  accurate 
in  time  and  space.  We  now  cover  second-order  accuracy  in  each. 

Second-order  time  accuracy  is  very  easily  achieved  by  simply  replacing 
the  two-point  backward  temporal  differencing  by  a  three-point  backward 
formula  in  time  which  for  equal  time  steps  can  be  written  as 


q 


t 


1.5 


n+1 

3 _ 


-  2  qn  +  0.5  qn_1 
__ 


(2.15) 


A  modified  formula  may  easily  be  obtained  for  varying  time  steps. 


Second -order  spatial  accuracy  is  a  iso  achieved  as  easily  as  time 
accuracy.  In  practice,  all  one  must  do  is  to  use  second-order  upwind 
discretizations  on  the  right  hand  side  of  the  equation  that  describes 
the  relaxation  algorithm  before  the  step  of  writing  the  Gauss-Seidel 
sweeps.  For  example,  the  x-derivative  term  may  be  discretized  as 


~  1  -5  q,  -  2  qj.)  +  °-5  q 


'  j  -2 


Ax 


(2.16) 


for  backward  difference  approximations  needed  with  positive  propagation. 
These  second-order  discretizations  may  have  to  be  suitably  modified  or 
limited  to  result  in  a  TVD  scheme  [Ref.  12].  Numerical  experiments  show 
that  even  without  limiting,  rapid  convergence  of  the  relaxation  iterations 
is  possible  for  the  problems  under  consideration. 


3.0  AX  I  SYMMETRIC  STREAM  FUNCTION  FORMULATION 


The  Nav ier-Stokes  equations  fit  very  well  into  the  framework  developed 
in  the  previous  section.  As  a  first  application,  we  consider  here  the  incomp 
ressible  flow  of  a  fluid  contained  in  a  spinning  shell.  For  simplicity  of 
treatment  here,  let  the  container  be  a  body  of  revolution  spinning  about  its 
axis.  If  z  were  the  axial  direction  and  r  and  0  the  polar  coordinates 
in  the  plane  perpendicular  to  the  axis,  the  governing  equations  may  be  exp¬ 
ressed  in  terms  of  an  axisymmetric  stream  function  \p  ,  a  vorticity  £  and 
a  circulation  y  (see  Ref.  3  from  which  the  governing  equations  are  adapted 
here)  in  the  r,  0,  z  coordinates  as 


2  *r 

-  -f  =  r  £ 

Ct  ♦  u  ?r  +  w  Cz  -  u  f 
Yt  +  u  Yr  +  w  Yz 


-  iv2;  *  4 


Re 


1 rn2  Tr 

■  R?IV  Y  -  -T 


Here, 

2  2  2 
V2  =*  — -  +  — „  ,  Re  =  0-JL_  s  Reynolds  number 

3r  9zZ  V 


Y  =  r  v  ,  u  =  (^2)/r  ,  w  =  -  (tp  ) /r 


(3.1a) 

(3.1b) 

(3.1c) 


(3. Id) 


where  ft  is  the  spin  rate  and  v  is  the  azimuthal  velocity. 

All  variables  above  are  non-dimensional.  The  non-dimensional ization 
used  is  given  below  (once  again  taken  from  Ref.  3). 


r  *=  R/a  ,  u  =  U/(Qa)  ,  v  =  V/(Qa)  ,  w  =  W/(fta)  ,  z  =  Z/a 
\p  =  ¥/(fla3)  ,  y  =  T/(na2)  ,  £  =  T/0  ,  t  -  «T 
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(3.2) 


The  dimensional  quantities  have  been  shown  above  as  capital  letters.  The 
normalizing  length  scale  a  is  usually  chosen  to  be  the  cylinder  radius. 

The  first  equation  ( E q .  3.1a)  is  an  elliptic  Poisson  equation  for  the 
stream  function.  The  other  two  ( E q .  3.1b  and  Eq.  3.1c)  are  parabolic  vor- 
ticity  transport  equations.  Without  the  viscous  terms,  these  latter  equa¬ 
tions  are  hyperbolic.  Note  that  these  two  equations  for  £  and  y  involve 
ther  terms  u  and  w  which  are  derivatives  of  the  stream-function,  as  only 
multipliers  of  derivatives  of  £  and  y  .  Also,  the  first  equation  involves 
no  derivatives  of  C  or  y  (in  fact  no  y  even  occurs  there).  Under  these 
circumstances,  it  is  conventional  procedure  to  first  update  vorticity  and 
circulation  from  the  n-th  time  level  to  level  n+1  .  Then,  the  updated 
vorticity  values  are  used  on  the  right  hand  side  of  the  equation  for 
stream  function  to  update  the  value  of  t^. 

3.1  Numerical  Method 

In  our  procedure,  along  with  this  implementation,  we  use  Gauss-Seidel 
iterations  for  all  three  equations  after  the  transport  terms  of  the  equations 
for  vorticity  and  circulation  are  approximated  using  one-sided  spatial 
derivatives  depending  on  the  sign  of  the  transport  velocity  multiplying 
them.  The  elliptic  part  (the  viscous  terms)  of  Eqs.  3.1b  and  3.1c  fit 
snugly  into  the  structure  of  Gauss-Seidel  Iterations  and  Eq.  3.1a  is  noth¬ 
ing  but  elliptic. 

We  now  write  down  the  discretizations  in  detail  for  Eqs.  3.1.  Let 
the  coordinates  r  and  z  be  indexed  with  the  subscripts  j  and  1 
respectively.  Thus,  j  and  1  also  denote  computational  coordinates 
with  grid  points  located  at  integer  values  of  j  and  1.  The  metrics 


are  then  given  by  3r/9j  ,  9z/91,  3  r/3j  ,  9  z/91  ,  etc.  The  following  trans 
formations  are  then  valid  for  any  quantity  f  : 


9f_  ,  9r 

9j  3j 


-  ( 


3  r /3r^9f^  ,  /  9  r  \  2 

^T9j;9T;M3j; 


9f_  .  3r 
3j  3j 


_  /9  f  /3  z.9zv3fv./3zv 

=  l^lT9T,9T)Arr) 


(3.3a) 


(3.3b) 


(3.*a) 


(3 .^b) 


Central  difference  formulae  for  any  quantity  f  are  given  by 

If  '  (fj*i  -  fj-.)/2  (3-5a) 

(92f/3j2)  =  fj+1  -  2  f  +  f  (3.5b) 


along  the  radial  computational  coordinate.  Similar  expressions  are  valid 
along  the  axial  coordinate.  Forward  and  backward  difference  formulae  for 
first  derivatives  are  given  by 


M 


(i 


5fj 


+ 


0.5 


(3.6) 


Velocities  u  and  w  are  computed  from  central 
Central  differences  are  used  for  all  derivative  terms 
derivative  terms  divided  by  the  Reynolds  number  (i.e. 


difference  formulae, 
in  Eq.  3.1a  and  a  I  I 
all  viscous  derivative 
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6 


Initial  conditions  are  given  by 


tf>(0,  r,z)  =  0,  Y(°.r,z)  =  (fl./fi)r2 


(3.10 


where  ft  is  the  final  cylinder  rotation  rate  (also  used  in  the  non-dimen¬ 
sional  izat  ion  process)  and  ft.  is  the  initial  rotation  rate  (  =  0  for 
spin -up  from  rest) . 
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4.0  PRIMITIVE  VARIABLE  FORMULATION 


The  primitive  variable  formulation  is  considered  here  in  r,  9,  z 
(cylindrical)  coordinates.  The  equations  and  method  considered  in  this 
section  will  be  used  in  Sections  6  and  7. 

4,1  The  Governing  Equations 

The  governing  equations  are  written  in  the  inertial  coordinate  system 
as  follows: 

Momentum  equations 

2 

9^u  +  u9ru  +  (v/r)9gU  +  w9zu  -v  /r  +  9rP 

9  v  +  u9  v  +  (v/r)9Qv  +  w9  v  +uv/r  +(l/r)9Qp 
t  r  u  z  t) 

9^w  +  u9rw  +  (v/r)9gW  +  w9^w  +  9zP 

Continuity  equation 

9  u  +  (l/r)9Qv  +  9  w  +  u/r  =  0  . 

r  W  z 

In  the  above. 


V2  = 

9^+0 

r)9r  + 

( 1/ r2) 9g  +  92 

(4.1c) 

and 

u  = 

rad ia  1 

velocity, 

r 

=  radial  coordinate, 

V  = 

azimuthal  velocity, 

9 

=  azimuthal  coordinate, 

w  * 

axial 

velocity, 

z 

=  axial  coordinate. 

The  same  non-d imens iona 1 i zat ions  employed  in  the  previous  section  has  been 
used  in  the  above  equations.  The  quantity  V  is  the  reciprocal  of  the 


=  v(V2u  -  (2/r2)9gV  -  u/r2) 

=  v(V2v  +  (2/r2)9gU  -  v/r2) 

=  v(V2w  ) 

(4.1a) 

(4.1b) 
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to1 


Reynolds  number. 


The  continuity  equation  has  no  time  derivative.  To  facilitate  const¬ 
ruction  of  an  implicit  numerical  method,  the  continuity  equation  is  modi¬ 
fied  (see  Refs.  6,  15,  16  for  a  similar  treatment  applied  to  approximately 
factored  schemes)  to  be 

3.p  +  6(9  u  +  (l/r)3Qv  +  9  w  +  u/r)  =  3  p*  .  (*».2) 

t  r  o  z  t 

The  quantities  6  and  p*  and  the  numerical  algorithm  will  be  chosen 
to  let 


9  u  +  (l/r)3„v  +  3  w  +  u/r  -*■  0 

r  t)  z 

Equations  4. la  and  U.2  may  be  combined  to  yield 


(**.3) 


A  Qr  +  (1/r)  B  Qq 


+  C  0 


where 

n  t  \  transpose  tJ  .  .  .  .  .  .  . 

Q  =  (u,  v,  w,  p)  r  ,  Hv  =  viscous  terms,  and  H.  =  inviscid 

terms  that  do  not  involve  derivatives,  and 


u 

0 

0 

r 

V 

0 

0 

*0 

w 

0 

0 

0~ 

0 

u 

0 

0 

,  B  = 

0 

V 

0 

1 

,  C  = 

0 

w 

0 

0 

0 

0 

u 

0 

0 

0 

V 

0 

0 

0 

w 

1 

e 

0 

0 

0 

0 

6 

0 

0 

0 

0 

6 

0 

(V5) 


4.2  The  SCM  Method 


The  Split  Coefficient  Matrix  (SCM)  method  of  discretization  is  applied 
to  the  hyperbolic  part  of  the  modified  Navier-Stokes  equations.  A  detailed 
description  of  the  SCM  method  can  be  foundin  Refs.  13  and  14.  The  method 
applied  to  the  equations  under  cons iderat ion  here  is  sufficiently  described 
now  for  completeness. 

Each  of  the  coefficient  matrices  A,  B,  C  can  be  written  as  (in 
the  following,  we  do  not  type  an  underscore  for  the  matrices  and  we 
delete  the  arrow  on  top  for  the  vectors) 


A  RA^ALA  ,  B  RBABLB  ’  C  RC^CLC  *  (**•&) 

Here,  the  subscripts  A,  B,  and  C  for  matrices  R,  A,  and  L  are  used  to 
denote  what  coefficient  matrix  the  latter  correspond  to.  The  matrix  of 
right  (or  column)  eigenvectors  of  a  coefficient  matrix  has  been  denoted 
by  R.  That  is,  each  column  of  R  is  a  right  eigenvector  of  the  corres¬ 
ponding  coefficient  matrix.  Similarly,  each  row  of  L  is  a  left  (or  row) 
eigenvector  of  the  coefficient  matrix.  The  diagonal  matrix  of  eigenvalues 
of  the  coefficient  matrix  has  been  denoted  by  A.  Thus,  for  example. 


A  ra  =  ra  aa  ’  la  a  =  aa  la  •  {k-7) 

In  Eqs.  4.6,  it  has  been  tacitly  assumed  that  the  left  and  right  eigen¬ 
vectors  have  been  suitably  normalized  such  that 

R  L  =  L  R  =  I  (4.8) 

where  I  is  the  identity  matrix.  When  the  coefficient  matrices  belong 


to  a  hyperbolic  system  of  equations,  the  left  and  right  eigenvector  matrices 
are  of  full  rank  (i.e.  within  each  of  these  matrices,  each  eigenvector  is 
linearly  independent  of  the  others). 

For  the  SCM  method,  the  eigenvalues  are  split  into  those  that  are  posi¬ 
tive  and  those  that  are  negative  (and  those  that  are  zero  can  be  in  either 
group) . 

A+  =  (A  +  |A|)/2  ,  A"  =  (A  -  i A ) ) /2 

Now,  the  split  coefficient  matrices  may  be  defined  to  be 

+ 

A±  ,  B±  ,  C±  =  RA jB>c  AA,B,C  LA,B,C 

Rewriting  each  coefficient  matrix  as  the  sum  of  its  positive  and  nega¬ 
tive  parts,  it  is  appropriate  to  use  backward  difference  approximations  for 
those  first  derivative  terms  that  are  multiplied  by  the  positive  part  of 
the  coefficient  matrix  and  vice-versa.  Thus, 

A  Qr  -  \  (irb  +  A.  Qrf  (A. 11) 

with  subscripts  b  and  f  denoting  backward  and  forward  differences. 
Equation  A . I  1  may  be  rewritten  as 

A  Qr  =  j(A+| A| )  Qfb  +  i(A-|A|)  Qrf  (A. 12) 

with  | A |  =  ]Aa|  L^.  Equation  A. 12  may  be  rewritten  in  turn  as 

A  Qr  =  A  (Qrb+Qrf)/2  +  |A|  (Qrb-Qrf)/2  (A. 13) 


(A. 9) 


(A. 10) 
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Qrb  =  (1.5  Qj  -  2  Qj_)  +  0.5  Qj.^/Or/ajJj 
Qrf  *  -0.5  Qj  -  2  Qj  +  ]  +  0.5  Qj+2)/(3r/3j)j 

Substituting  this  into  Equation  4.14,  we  obtain 


A  ar  *  (57737)  Aj<V,-Qj-r<Qj.2-Qi-2)/'' 


+  (3773T) 1  sj  I ( Qj  -2 ■l|Qj  - 1 +6Q i  ■1| ti ^^+2’ /l* 


(4.15) 


(4.16) 


Matrix  A: 


The  first  term  on  the  right  hand  side  is  a  central  difference  approxi¬ 
mation  (albeit  an  unusual  one)  to  AQ^ .  The  second  term  is  a  fourth 
difference  diffusion  term  that  arises  naturally  in  the  SCM  method  and 
helps  squash  high  frequency  oscillations. 


The  various  eigenvalues  and  matrices  are  now  defined.  For  coefficient 


(4.17a) 


AIA  =  u  «  A2a  =  u  >  A™  =  (u+(u2+4B)  1/2)/2  ,  X.  =  (u-(u2+46)1/2)/2 


(uN-46) 


(4.17b) 


(u2+46) 1/2 


~a4a 

(u2+46) 


(u2+46),/2 
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The  eigenvalues  and  eigenvector  matrices  for  the  coefficient  matrices  B  and 
C  are  similarly  derived. 


A. 3  Relaxation  of  the  Discretized  Equations 


The  SCM  method  given  above  is  used  to  discretize  the  inviscid  terms. 
Central  difference  approximations  are  used  for  the  viscous  terms.  Three- 
point  or  two-point  backward  discretization  of  the  time-derivative  terms 
leads  to  second-order  or  first-order  time  accuracy  as  desired.  The  time- 
discretized  implicit  version  of  Eq.  k.k  can  then  be  written  as 


(l+4>/2)( 


(TV) 


,Qn-Qn-] 


(<p/2)  (*■  g  )  +  (AQp+  ( 1  /r ) BQq+CQz+H . ) 


(A .20) 


+  (0,0,0, 3tp*) 


transpose 


with  (}>  =  0  for  first-order  time  accuracy  and  4>  =  1  for  second-order 
t ime  accuracy. 

Writing  q’  to  be  the  i -th  iterate  approximation  to  Qn+' ,  we 
can  write  the  Gauss-Seidel  method  for  Eq.  ^.20,  coupled  with  appropriate 
SCM  (upwind)  and  central  difference  approximations,  to  be 

3H  3H 

!>' '*3T75TlAl4  imzW+jmrM- -‘sr-f'sf-r1" (^i-qi',) 

J,K,I  J  »  K  ,  I 


0^/2)jr(qi",-Qn)-  +  <AQr  *  X  ♦ 


('*.21) 


+  (H. )'*'■'  -  +  (0,0,0,p*)tran 


spose 


To  implement  the  Gauss-Seidel  sweep,  all  terms  with  the  superscript 


.  *•  V.v  -V 


.A  .  «__*V  -  .  V.  ,  . 


are  evaluated  using  all  available  i-th  level  values  along  with 
i-1  level  values  for  those  quantities  for  which  i-th  level  values  have  not 
been  computed  in  the  sweep.  Each  cycle  of  sweeps  comprises  a  forward  sweep 
(increasing  values  of  indices  j,  k  and  1)  and  a  backward  sweep  (decreasing 
values  of  indices,  starting  from  the  maximum  values  of  these  indices). 

The  pseudo  pressure  transient  term  p*  in  Eq.  4.21  is  taken  to  be 
P*  -  (1+  |)(pM-pn)/(At)  -  (4>/2)  (pn-pn"1)/(At)  (4.22) 

Thus,  when  the  relaxation  iterations  converge  between  two  time  steps,  the 
two  transient  terms  for  pressure,  namely  (p'-pn)  and  (p'  ^-pn)  will 
cancel  each  other  out,  leading  to  exact  satisfaction  of  the  continuity 
equation.  However,  perfect  convergence  of  the  relaxation  iterations  is  not 
sought  for  every  time  step;  approximate  convergence  of  the  subiterations 
will  almost  always  suffice. 

4.4  Advantages  of  the  Unfactored  Scheme 

Approximately  factorized  schemes  incur  an  error  due  to  the  factoriza¬ 
tion.  Also,  for  the  imcompress ible  Navier-Stokes  equations,  when  the  conti¬ 
nuity  equation  in  modified  as  in  Eq.  4.2,  for  large  values  of  B,  the  momentum 
equations  are  contaminated  by  this  factorization  error.  Unfactored  schemes 
avoid  this  error.  Thus,  larger  values  of  B  can  be  chosen.  When  B  ~  1 / (At) 
the  exact  continuity  equation  (Eq.  4.1b)  is  satisfied  to  first-order  accu¬ 
racy.  And  larger  values  of  3  lead  to  a  greater  degree  of  satisfaction  of 
the  continuity  equation.  The  other  advantage  of  unfactored  implicit  schemes 
is  that,  as  a  steady  state  is  approached,  larger  and  larger  values  of  the 
time  step  may  be  taken. 
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5.0  ROTATING  FRAMES  OF  REFERENCE 


Thus  far,  we  have  been  considering  an  inertial  reference  frame  and  the 
dependent  variables  were  the  inertial  velocity  components.  It  is  often  very 
convenient  to  work  with  a  frame  of  reference  attached  to  the  spinning  and 
precessing  shell  and  also  compute  the  velocity  components  in  the  rotating 
frame  of  reference.  The  governing  equations  are  given  below  for  a  coor¬ 
dinate  frame  attached  to  a  rotating  shell.  The  treatment  follows  that  of 
Greenspan  [Ref.  17] . 

Let  L,  U  U  characterize  reference  length,  time  and  velocity, 
respectively,  of  a  particular  motion.  The  position  vector  7 ,  time  t, 
velocity  vector  q,  rotation  vector  ,  and  pressure  can  then  be  scaled 
respectively  by  L,  fl  *  ,  U,  ft  and  pftUL.  The  vectors  are  in  a  Cartesian 
frame  of  reference  fixed  to  the  spinning  and  precessing  shell.  Then  for 

(t)  =  il  (  k  +  e  J(t)  )  (5.0 

the  governing  equations  may  be  written  in  non-dimensional  form  as 

Momentum  equations: 

3tq  +  eq.^q  +  2 ( k  +  c^(t))  x  q  =  -^p  +  r  x  (d/dt)^  +  EV2q 
Continuity  equation: 

$.q  ■  0 

where  the  Ekman  number  is  given  by 
E  -  v/(«L2)  (5.3a) 

and  the  Rossby  number  is 


(5.2a) 

(5.2b) 
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IF.  LFJ  V"  F'.'  i  J 


A 


i  ■  r  ■  y  *  u  ■  »  ■  **  pi  i  «  u; 


rTr  ■  f  v  »7  rrv 


e  -  U/(flL)  (5.3b) 

The  reduced  pressure  p  appearing  above  inc 1 udes  along  with  the  static 
pressure,  any  body  forces  (assumed  to  be  conservative),  and  the  centri¬ 
fugal  acceleration. 


p  =  static  pressure  +  conservative  body  force 
-(l/2)p(ft  x  r ) . (fi  x  r) 


(5.4) 


The  z  axis  and  consequently  the  k  unit  vector  is  assumed  to  be  the 
axis  of  spin  of  the  container.  The  precession  £(t)  is  superposed  on  this 
basic  spin. 


The  governing  equations  5.1  and  5.2  have  been  written  in  a  Cartesian 
coordinate  system  attached  to  the  spinning  shell.  The  z  coordinate  lies 
along  the  spin  axis.  The  following  remarks  are  useful  to  transform  these 
equations  into  a  set  of  equations  written  for  a  cylindrical  coordinate  system 
attached  to  the  rotating  shell  and  for  cylindrical  dependent  variables  in 
that  system.  The  z  axis  is  identical  in  both  Cartesian  and  cylindrical 
systems.  The  other  axes  are  linked  by 


x  =  r  cosQ  ,  y  =  r  sin8 


(5.5) 


The  velocity  components  are  linked  by 


Ucartesian  = 
v 

cartesian  ** 
w 

cartesian  = 


u  cos0  -  v  sin9 
u  sin9  +  v  cos9 

w 
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(5.6) 


where  the  cylindrical  velocity  components  have  been  denoted  as  usual  by 
u,  v  and  w. 

The  metrics  transform  according  to 


X  = 

cos9  , 

X- 

=  -r  sin0 

r 

9 

(5.7) 

yr  = 

s  in0  , 

y0 

=  r  cos0 

and 

r  * 

COS0  , 

0 

=  -  (sin9)/r 

X 

X 

(5.8) 

r  = 

s  in0  , 

0 

=  (cos0)/r 

y 

y 

There  are  three  scalar 

equations  identifiable  in  the  vector  equation 

5.2a, 

one  each  for  each 

Cartesian  velocity  component.  To  go  from  the  equa 

t  ions 

in  5.2a  and  5.2b 

to  a 

cylindrical  variable  formulation  similar  to 

Eq.  k 

.  la  and  Eq.  k . lb, 

we  can  take  the  following  steps.  First,  in  each  of 

the  equations  in  5.2,  transform  the  independent  variables  using 


3/3x  = 

cos0  3/3r 

(l/r)sin0  3/39 

3/3y  = 

sin0  3/3r 

+  (!/r)cos0  3/30 

(5.9) 

32/3x2  - 

2  2 

►  3Z/3y  = 

32/3r2  +  (l/r)3/3r  +  (l/r2)32/302 

Now,  we  have  four  equations  for  the  Cartesian  velocity  components 
but  in  the  cylindrical  independent  variables.  Into  each  of  these  equa¬ 
tions,  substitute  Eqs.  5.6  to  obtain  four  equations  involving  the  cylind¬ 
rical  velocity  components.  The  continuity  equation  and  the  equation  for 
Cartesian  w  component  of  velocity  will  now  be  in  the  form  of  Eqs.  *».!. 


However,  the  other  two  equations  will  each  involve  linear  combi nat ions 
of  u  and  v.  To  obtain  mutually  exclusive  prediction  equations  for 


u  and  v,  recombine  these  two  equations  as  follows:  to  obtain  an  equation 


for  u,  multiply  the  equation  for  u  .  by  cos9,  the  equation  for 

ca  ries lan 

v  .  by  sin9  ,  and  add;  to  obtain  an  equation  for  v,  multiply  the 

cartes lan 

equation  for  u  „  .  by  (-sin9),  the  equation  for  v  ...  by  (cos9) 

cartesian  1  *  ^  cartesian  1 

and  add . 


This  procedure  has  been  used  for  the  sets  of  governing  equations  given 
in  the  next  two  sections.  The  difference  between  those  sets  and  Eqs.  4.1 
can  be  summarized  as  follows:  all  convective  terms  are  now  scaled  by  the 
Rossby  number;  new  source  terms  corresponding  to  the  Coriolis  acceleration 
2 ( k  +  ej(t))  x  q  and  corresponding  to  the  acceleration  of  the  rotating 
frame  given  by  "r  x  d6/dt  ,  are  added.  Thus,  it  is  an  easy  matter  to 
adapt  the  methodology  developed  for  the  primitive  variable  formulation  in 
the  last  section  for  the  inertial  coordinate  system  to  rotating  frames  of 
reference.  The  source  terms  have  no  spatial  derivatives  in  them  and  are 
added  in  a  very  straight-forward  manner.  For  computing  efficiency,  they 
must  however  be  treated  implicitly  in  the  relaxation  iterations.  That  is, 
the  Jacobian  matrix  of  these  source  terms  must  be  included  along  with  the 
other  terms  in  the  diagonal  matrix  of  the  left  hand  side  of  the  operator  in 
the  implicit  method. 
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6.0  CONING  -  SPIN  AND  PRECESSION 


We  now  apply  the  governing  equations  developed  in  the  last  section  for 
rotating  frames  of  reference  to  the  particular  case  of  simple  coning  which 
is  a  combination  of  spin  and  precession.  Let  the  cylindrical  shell  be  spin 
ning  with  spin  rate  ft  and  precessing  with  precession  rate  w.  We  define 

Tnu  -  a)  /  ft  (6.1) 

from  which  we  can  define 

a  «  1/  0  +  Tnu)  (6.2) 

Let  the  axis  be  precessing  with  a  yaw  angle  a  .  In  a  slight  departure 
from  the  normalization  procedure  of  the  previous  section,  we  define  first 
the  Rossby  number  to  be 

e  =  a  a  (6.3) 

from  which  we  can  derive  the  normalizing  velocity  U  to  then  be 

U  =  e  ft  L  (6.4) 

where  L  is  the  usual  length  scale  (typically  taken  to  be  the  radius  of 
the  cylinder).  All  we  have  done  here  is  to  pick  the  Rossby  number  first 
and  obtain  the  non-d imens iona 1 iz ing  velocity  from  it,  rather  than  the  other 
way  around. 

The  perturbation  £(t)  is  then  given  by  [Ref.  17] 
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6.5 


$  =  -Tnu  (cos (at)  i  +  sin (at)  j  )  ( 

The  source  terms  in  Eq.  5.2a  are  thus  given  by 


2 (k  +  E^(t) )  x  q  = 


and 


x 

-O 

II 

-  (v  .  + 

cartesian 

eT 

NU 

+  (u  .  .  + 

cartesian 

eT 

^NU 

+  eT....(u  .  .  sinat  +  v  .  .  cosat)  k 

NU  cartesian  cartesian 


z  cos  at  i 


x  d^/dt 


*  aTNU  +z  sin  at  j 


-(x  cos  at  +  y  sinat)k 


(6.6) 


(6.7) 


In  the  cylindrical  rotating  frame  of  reference,  the  source  terms  that 
must  be  added  to  the  left  hand  side  can  be  computed  from  the  above. 


Source  term  for  radial  momentum  equation 

-  -2v  -  2eTN|jSin(at-0)  -  zaT^cos (at-9)  (6.8) 


Source  term  for  azimuthal  momentum  equation 

*  +2u  +  2eTNUcos (at-6)  -  zaTN(Jsin(at-0)  (6.9) 


Source  term  for  axial  momentum  equation 

*  2eT.M1[u  sin(at+0)  +  v  cos(at+9)]  +  aT^-.fx  cosat  +  y  sinat] 
NU  NU 


The  boundary  conditions  along  all  walls  are  given  to  be 
u  ■  v  “  w  *  0 


in  the  rotating  frame. 


(6 


(6.1 
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7.0  SPINNING  SHELL  WITH  PRECESSING  LID 


In  this  section,  we  consider  a  spinning  cylindrical  shell  with  only  the 
lid  precessing  at  one  end  of  the  cylinder.  Once  again,  a,  ft,  w  are  defined 
to  be  the  precession  angle,  spin  rate  and  precession  rate,  respectively  (as 
before).  However,  for  this  case,  we  define 


a  =  1  -  T, 


(7.D 


and  we  define  the  Rossby  number  to  be 


c  =  sin (a) 


(7.2) 


The  rotating  frame  of  reference  does  not  precess.  Thus, 


(7.3) 


.  ..  .  ....  .  .  This  contributes  -2 v 

and  the  source  term  is  simply  given  by  2k  x  q 

to  the  radial  momentum  equation  and  2u  to  the  azimuthal  momentum  equation. 

The  boundary  conditions  along  the  walls  are  given  as  follows: 

Along  the  walls, 


u  =  v  ”  0 


(7.4) 


Along  the  cylinder  wall  and  lower  end  wall. 


w  =  0 


(7.5) 


Along  the  upper  end  wall, 

w  ■  r  a  sin(ot+0^) 


(7.6) 


ire  8^  is  a  given  phase  angle. 
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8.0  CONCLUDING  REMARKS 


A  new  formulation  is  presented  for  constructing  finite  difference  algo 
rithms  for  the  incompressible  Navier-Stokes  equations.  The  new  method  is 
based  on  constructing  relaxation  procedures  for  unfactored  implicit  upwind 
schemes  for  the  hyperbolic  part  of  the  system  of  governing  equations  taken 
together  with  central  differencing  for  the  viscous  terms.  The  new  approach 
has  been  applied  to  the  axisymmetric  stream  function  formulation  as  well  as 
for  the  primitive  variable  equations  in  three  spatial  dimensions.  The 
formulation  has  been  extended  for  rotating  frames  of  reference.  The  two 
cases  considered  here  are  1)  coning  motion,  and  2)  spinning  shell  with 
precessing  lid.  The  computer  programs  developed  in  this  study  have  been 
installed  at  the  Ballistic  Research  Laboratory,  Aberdeen  Proving  Ground, 
and  are  being  used  to  obtain  results  by  BRL  personnel  for  many  cases  of 
interest.  Those  results  will  be  incorporated  into  a  BRL  Technical  Report 


at  a  future  date. 
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